COVER SHEET 

NOTE: This coversheet is intended for you to list your article title and author(s) name only 
— this page will not appear on the CD-ROM. 


Paper Number: 1500 ( replace with your paper number ) 

Title: An enriched shell element for delamination simulation in composite 

laminates 

Authors: Mark McElroy 


(FIRST PAGE OF ARTICLE - align this to the top of page - leave space blank above ABSTRACT) 


ABSTRACT 

A formulation is presented for an enriched shell finite element capable of 
delamination simulation in composite laminates. The element uses an adaptive splitting 
approach for damage characterization that allows for straightforward low-fidelity model 
creation and a numerically efficient solution. The Floating Node Method is used in 
conjunction with the Virtual Crack Closure Technique to predict delamination growth 
and represent it discretely at an arbitrary ply interface. The enriched element is verified 
for Mode I delamination simulation using numerical benchmark data. After determining 
important mesh configuration guidelines for the vicinity of the delamination front in the 
model, a good correlation was found between the enriched shell element model results 
and the benchmark data set. 


Mark McElroy, NASA Langley Research Center, 2 W. Reid St, Mail Stop 188E, Hampton, VA 
23681 USA. 



INTRODUCTION 


Progressive damage simulation in composite laminates using finite elements is often 
performed using models consisting of solid elements with a mesh refinement of at least 
one element per ply. This level of high-fidelity is necessary with current simulation 
techniques to be able to predict and represent physically complex three-dimensional 
damage processes through the thickness of a laminate. Excessive computational demand 
and required user expertise for complex high-fidelity simulations, however, often limit 
the use of these tools in practice by driving up the time and cost associated with a given 
analysis. As the continuous advancement of computational technology allows for 
increasingly complex models and shorter run times, the cost of necessary user training 
and the required time for model creation and verification still may remain high. 

Model fidelity is typically high for two reasons. The first reason is the need for 
detailed stress and strain information to be available through each ply for a criterion to 
use for prediction of initiation and propagation of a laminate damage process such as 
matrix cracking or delamination. The second reason for high model fidelity is simply to 
have a mesh capable of representation of a three-dimensional damage pattern as dictated 
by the damage growth criterion. The price paid for high-fidelity models such as this, in 
general, is that they are very inefficient as a large amount of information is determined 
throughout the entire model when this level of detail is only needed at locations where 
a damage process is occurring. 

Shell elements have long been used by industry for analyses that are less complex 
than laminate damage simulation and have proven to be a cost effective analysis tool. 
Shell element models are computationally efficient, and due to a simplified mesh when 
compared to a high-fidelity solid element approach, the models are inherently less 
expensive for an analyst to create and use. Use of shell element models for laminate 
damage simulation, however, introduces a number of challenges, the most obvious of 
which may be prediction and characterization of out of plane transverse matrix cracks 
and delaminations at multiple interfaces. A fundamental prerequisite for modeling a 
complex damage process such as this is the use of shells to simulate simple delamination 
problems that do not include matrix crack interaction. This capability should be 
demonstrated before further work towards use of shells for a more complex three- 
dimensional simulation is undertaken. 

Layer-wise stacking is a common technique for using shell elements for laminate 
damage simulation [1,2]. In such an approach, independent layers of shell elements are 
defined and connected across ply interfaces in a laminate. If delamination at a specific 
ply interface is of interest, a mesh may consist of two layers of shell elements, one on 
either side of the ply interface location. The two layers can be connected by rigid links 
and delamination growth can be simulated using the Virtual Crack Closure Technique 
(VCCT) [3,4], Wang et al. [5] were among the first to study the use of VCCT for 
delamination growth in this manner and presented formulas for calculation of energy 
release rates using shell elements. This work, along with a related study later [6], 
indicate that mixed-mode energy release rate can be calculated using shell elements, 
however, continued mesh refinement may not result in converged values in some cases. 
Similarly, in bimaterial crack tip interfaces, energy release rate component divergence 
is observed during mesh refinement when using VCCT with any element type. Previous 
work indicated that this non-convergence may be avoided by sizing the mesh above a 



minimum value based on the ratio of crack growth length (i.e., mesh size) to element 
thickness, Aa/h [7], 

An inherent disadvantage of the two-layer stacked shell models is that the 
delamination interface must be predefined, thereby limiting the model’s utility as a 
general predictive tool. Furthermore, by restricting the delamination to a single ply 
interface location, multiple delaminations, or delamination and matrix crack interaction, 
is not possible. One means to avoid this is to model many or all interfaces in the mesh 
where there is an independent layer of shells for every ply. This type of approach, 
however, begins to share some of the same drawbacks mentioned for high-fidelity solid 
element models. Ideally, in terms of computational efficiency and ease of use, a 
laminate would be modeled as a single layer of shell elements where a delamination 
could form and propagate at an arbitrary location in the layup. This approach has been 
enabled for through thickness cracks by numerical simulation techniques such as the 
Extended Finite Element Method (XFEM) [8], the Phantom Node Method (PNM) [9], 
and the Floating Node Method (FNM) [10]. The same concept can be applied to 
delamination in shell elements, but the arbitrary discontinuity is constrained to be 
orthogonal to the shell surface normal. The fidelity required for damage growth in a 
model such as this can be thought of as adaptive in that the model is defined initially in 
low-fidelity, one element thick, and remains in this state everywhere except locally 
where delamination occurs and multiple layers are required. Simulation models based 
on shell elements that allow for an arbitrary delamination location and adaptive fidelity 
through the thickness have been proposed and studied only recently [11-13]. 

The goal of this paper is to present the formulation and initial verification of a novel 
enriched shell element capable of delamination simulation in composite laminates. The 
element enrichment allows for adaptive mesh fidelity through the thickness where a 
single element splits into two elements only as required locally to characterize an 
evolving delamination process. At any location in the model where damage does not 
occur, the original discretization remains unaffected and a single shell element is used 
to represent the entire laminate thickness as defined in the beginning of the analysis. 
This capability is achieved using the Floating Node Method [10] to discretely represent 
a delamination in the mesh and YCCT to predict its growth. The element is coded as a 
user-subroutine in Abaqus 6.14/Standard® [14]. Details concerning the use of VCCT in 
conjunction with shell elements, mesh dependency of the solution (i.e., Aa/h ratio), and 
delamination fronts misaligned with the mesh are investigated. This study serves as the 
first milestone of the enriched shell element’s development process where the ultimate 
goal is to include the capability to simulate mixed-mode delaminations and transverse 
matrix crack-delamination interactions. 


ELEMENT FORMULATION 


Mindlin Shell Element 

The enriched finite element is based on a linear four node, six degrees of freedom 
(DOF) per node, Mindlin shell as presented by Onate [15]. The element has the ability 
to undergo transverse shear deformation and performs well for both thick and thin shells. 
Transverse shear correction factors are determined layer- wise as proposed by [16]. 



Shear locking is eliminated by imposing an assumed linear distribution of transverse 
shear strain in the element as described by Onate [15]. Full four point Gaussian 
integration is used to calculate the stiffness matrix. Classical laminate theory is used to 
determine orthotropic bulk material properties for the element based on a layup of 
differing ply orientations. 


Implementation of the Floating Node Method (FNM) 


The FNM works by embedding extra unused nodes as “floating DOF” into an 
element stiffness matrix as rows and columns of zeroes. If and when an element is split 
(see VCCT criteria in the next section), the floating DOF are activated and used to define 
what is effectively an additional element but still within the original element definition. 
Numerically, the FNM functions as follows. The stiffness matrix, fG e) , without any type 
of discontinuity in the displacement field, is determined by integration across the 
element area in the x- and y-directions and the element thickness, h, in the ^-direction 
and is given by 
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where H, a function of four linear shape functions, is the strain-displacement matrix, C 
is the constitutive material matrix, and b is the edge dimension of a square element. In 
a laminate shell formulation, C may be defined as a function of submatrices A, B, D, 
and G known from classical laminate theory that correspond to membrane, membrane- 
bending coupling, bending, and transverse shear behavior respectively. The A, B,D, and 
G matrices for a laminate are determined by 
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where C p and C s correspond to planar and transverse shear constitutive material 
properties, respectively. Equation (1) for a laminate shell can be re-written as 
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When a discontinuity is introduced at a single uniform z-coordinate, as in the case 
of a delamination, the element material is split into two regions, Qa and Qb, where 
region Qb corresponds to the DOF of the floating nodes. The stiffness matrix for a split 
element is given by 
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where K n ^ and /Q are determined using equation (3), but with their respective 
material integration limits split along the z-axis according to the discontinuity location. 
A diagram of an element in a damaged and undamaged state with associated stiffness 
matrices is shown in Figure 1. While the FNM has been used before in combination 
with VCCT [17], it has yet to be implemented in shell elements for in-plane 
delamination damage modeling as described. 

Currently, the FNM is implemented in the enriched element such that one 
delamination is possible, however, multiple delaminations could be handled with the 
same approach. As such, the element is defined with eight nodes total, four “real nodes” 
used in the definition of the initial undamaged region and four “floating nodes” used in 
the event that the element is split. There are 6 DOF per node resulting in a total of 48 
DOF per element. Computational efficiency is achieved by only performing stiffness 
integration across “activated” regions thereby numerically treating undamaged elements 
as if they were non-enriched conventional shell elements. This concept is illustrated by 
the stiffness matrix formulas shown in Figure 1. Unused floating nodes, as rows and 
columns of zeros, are removed from the global stiffness matrix by the Abaqus® solver 
reducing total DOF and effectively rendering undamaged elements back to a non- 
enriched state for the solution of equilibrium equations at each load increment. 
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Figure 1. Shell element damage states. 


VCCT and Damage Progression 

The enriched shell element can be used in a mesh to represent four different damage 
states. These are illustrated in Figure 2a and consist of 1) undamaged element, 2) split 
element ahead of the crack tip with all nodes tied, 3) split element at the crack tip with 
one or more nodes free, and 4) split element in the crack wake with all nodes free. At 



the end of every solution increment, VCCT is used on each delamination front tied node 
to determine if the tie should remain in place or be released resulting in a change in the 
adjacent element damage states (same as [18]). The tied nodes are constrained by 
inserting stiff springs internal to an element stiffness matrix connecting real and floating 
DOF (spring stiffness = 10 6 N/mm). 



crack extension length, Aa = b 
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b) Plan view of staggered delamination front with possible growth paths that the 
simulation considers at a corner crack tip node ij. 

□ = undamaged element 

I = split element ahead of crack tip w/all nodes tied 

□ = split element at crack tip w/one or more nodes tied 

□ = split element in wake of crack tip w/all nodes free 
— = delamination front 

All nodes have 6 global DOF u, v, w, 6 x , 0 y , 6 Z associated with nodal 
forces F x , F y , F z , M x , M y , M z respectively 

Figure 2. Dela mi nation front mesh description. 
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(i, j) = node designation 



In general, the use of VCCT is relatively straightforward when the mesh is aligned 
with the delamination front. However, when the delamination front and mesh are 
misaligned, the application of VCCT may not be intuitive. Methods have been proposed 
to implement VCCT such that it can be used accurately to simulate delamination growth 
independent of mesh alignment. One such approach involves determination of a growth 
direction vector which is used to calculate the required crack growth area and the 
location where the crack opening displacements are obtained [19,20]. Alternatively, 
Orifici et al. [18] presented an approach where correction factors are applied at 
delamination front nodes depending on the connectivity of the neighboring nodes. The 
approach used in the current study is less elaborate than any of these methods. 

Figure 2b can be seen as an example of a small section of a staggered delamination 
front where globally neither the shape of the front nor the direction of propagation is 
aligned with the mesh. VCCT is applied at the comer node (ij) to calculate the total 
energy release rate, Gt- The application of VCCT in the enriched shell element is based 
on an assumption that at a given tied crack front node, the delamination can only 
advance in orthogonal directions along mesh lines and therefore always obtains wake 
displacements at nodal locations. At the end of an increment, the energy release rate at 
a tied node is calculated in four directions parallel to the mesh: +x, -x, +y, and -y as 
shown in Figure 2b. Using delamination growth in the positive x-direction as an 
example, Mode I, Mode II, and total energy release rates at node (ij) are determined 
using formulas from [5] given as 
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where F, M, w, u, and 9 are nodal force, moment, z-displacement, x-displacement, and 
rotation, respectively. Nodal designations with a prime superscript refer to the “upper” 
set of elements (i.e., floating nodes) and crack extension area, AA (e> = b 2 , is the area of 
an element for a square regular mesh (in the current form of the VCCT implementation 
in the element enrichment, the mesh is required to consist of regular square elements). 
The critical energy release rate, G c , is obtained using the Benzeggagh-Kenane criterion 
[ 21 ] 


G? X) = G IC + (G //c -C /c )(g / ( / +x V4 +X) ) ??BK ( 8 ) 

where Gi c and Gu c are fracture toughness properties of the material and t|bk is a curve fit 
parameter determined from testing. 

Gt and G c are calculated in the other three propagation directions, -x, +y, and —y, by 
changing the wake displacement location, ( i-l,j ), in equations (5) and (6) to ( i+lj ), (i,j- 
1), and ( i,j+l ), respectively. The maximum Gt, out of the four possible propagation 
directions, is then compared to G c and if Gt > G c , the tie is released. Even though the 



enriched shell considers four possible growth directions, those that do not have any 
physical meaning (i.e., such as “reverse” growth in the -x-direction in Figure 2a, etc.), 
produce Gt = 0. If a delamination front node has more than one non-zero Gt, such as 
the comer node (ij) in Figure 2b, only the maximum Gt has any bearing on whether or 
not that tie is released. 

Mesh Configuration Near the Delamination Front 

It was found that when modeling a double cantilever beam (DCB) specimen the 
constrained DOF in the ties at the delamination front and the number of split (but still 
tied) elements ahead of the delamination front both affect the element performance 
significantly. Therefore, three cases were considered as summarized in Table I. In Case 
3, the number of split elements is based on a distance ahead of the delamination front 
and therefore varies between mesh sizes. The number of split elements ahead of the 
delamination front is maintained as the front advances in the simulations. 


TABLE I. SUMMARY OF THREE CASES CONSIDERED VARYING DELAMINATION FRONT 
TIED DOF AND NUMBER OF ELEMENTS SPLIT BUT TIED AHEAD OF THE DELAMINATION. 


Case 

DOF tied 

Number of elements 
split ahead of crack tip 

1 

lly V, Wy Oxy Qyy @z 

1 

2 

U, V, w 

1 



7, mesh =1.0 mm 

3 

U, V, W 

4, mesh = 2.5 mm 



2, mesh = 5.0 mm 


Due to the symmetric nature of the DCB specimen configuration and loading, each 
bending arm can be thought of as a cantilever beam supported at the crack tip with a 
point load at the end. The boundary condition at the crack tip is somewhat complex as 
the material supporting the bending arms ahead of the crack tip can deform in the axial 
direction and effectively allow a small amount of rotation. This root rotation would be 
captured in a high-fidelity model where the specimen is meshed through the thickness 
with multiple elements, however, this level of mesh fidelity is not present in a beam or 
shell model. Cotterell et al. proposed a method to simulate a small amount of root 
rotation in a simple beam model by increasing total crack length (i.e. increasing 
flexibility) [22], but it is not immediately obvious if this type of approach is viable in a 
general three-dimensional model. The three Cases described in Table I are investigated 
in an attempt to determine an accurate means to include the effects of root rotation in 
the structural stiffness of the specimen. 

Case 1 has rotations constrained in the ties at the delamination front and therefore 
does not account for any root rotation. Case 2 is considered as means to ‘relax’ the 
boundary condition present in Case 1 and allow some form of root rotation. Case 1 and 
Case 2 are compared in Figure 3 using two beam models representative of the upper 
DCB bending arm in each case. By summing forces and moments and comparing 
analytical beam deflection formulae, it is shown that in Case 2, the tie forces located at 
the crack tip and beam deflections to be used in VCCT calculations are a function of 
mesh size, b. 



a) rotations fixed at crack tip (Case 1) 


b) rotations free at crack tip (Case 2) 
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Figure 3. VCCT force and deflection dependency on constrained DOF at crack tip for 

Case 1 and Case 2. 


Case 3 is considered as a means to correct the mesh dependency in Case 2 while 
still providing a ‘relaxed’ boundary condition when compared to Case 1 that effectively 
allows some root rotation. A diagram of the Case 3 mesh is shown in Figure 4 where 
multiple elements are split ahead of the delamination front along the axis a up to a 
distance of O* max* 
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Figure 4. Diagram of mesh for Case 3: rotations are unconstrained in split 
element ties and multiple elements ahead of the delamination are split and tied 

together. 





ANALYSIS 


Results from a numerical Double Cantilever Beam (DCB) benchmark study [23] 
were used to verify the enriched shell element’s capability to simulate Mode I 
delamination. All analyses were static with linear ramping loads and ran in Abaqus 
6.14/Standard® [14]. The enriched shell element was coded in a user defined element 
subroutine (UEL). Models were created using the enriched shell element with three 
mesh sizes for Cases 1-3. The mesh sizes considered, 1.0 mm, 2.5 mm, and 5.0 mm, 
correspond to Aa/h ratios of 0.66, 1.33, and 3.33, respectively, where Aa is crack 
extension length and h is the thickness of one region of a split element. A schematic of 
the DCB test specimen and an overview of the finite element model are shown in Figure 
5. Because the model consists of planar coincident shell elements, it is shown in a 
deformed state to better illustrate its features. Material properties for the specimen are 
given in Table II [23]. Only Mode I verification is presented in this paper, though work 
is ongoing regarding Mode II and mixed-mode delaminations. 


TABLE II. MATERIAL PROPERTIES FOR DCB SPECIMEN: T300/1076 [23], 


Parameter 

Value 

Units 

En 

139.4 

GPa 

E22 

10.16 

GPa 

Ess 

10.16 

GPa 

Gl2 

4.6 

GPa 

Gis 

4.6 

GPa 

G23 

3.54 

GPa 

V12 

0.3 

- 

Vl3 

0.3 

- 

V23 

0.436 

- 

Gic 

0.17 

kJ/m 2 

Guc 

0.494 

kJ/m 2 

V]BK 

1.62 

- 




T300/1076: [0 24 ] Boundary condition 
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a) Specimen configuration 


Figure 5. DCB specimen. 


RESULTS & DISCUSSION 


Force displacement correlations are shown in Figure 6 for Cases 1-3. In all cases, 
the “sawtooth” behavior seen in the force displacement curves can be explained by the 
sudden release of a tied crack tip node as the delamination advances one element length. 
This behavior is seen to diminish with mesh refinement. In Case 1, based on the close 
match for the critical force magnitude seen in Figure 6a, the energy release rate appears 
to be predicted well but the stiffness of the model prior to damage is over predicted by 
approximately 24%. Also of note in the Case 1 results is that no mesh dependency is 
seen concerning stiffness or energy release rate prediction. Both of these observations 
are in accordance with the beam analysis previously shown in Figure 3 a. 

Case 2 does not constrain rotations in the ties and splits one element ahead of the 
delamination. It is clear from Figure 6b that while it is possible to match the stiffness 
well with this approach, both the stiffness and energy release rate are mesh dependent. 
These observations also are in accordance with the beam analysis shown previously in 
Figure 3b. 

Results from the Case 3 configuration, where tie rotations are unconstrained and 
more than one element ahead of the delamination is split, are shown in Figure 6c. Case 
3 was found to produce accurate stiffness and energy release rate results and remove the 
mesh dependency seen in Case 2. The applied force at a given displacement for Case 3 
generally is under predicted by approximately 5% for all mesh sizes, although at the 
critical value where damage initiates, the error is <1%. 



a) Case 1: Fixed rotations at crack tip 
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c) Case 3: Free rotations at crack tip - mesh split 
multiple elements ahead of delamination front 


Figure 6. DCB: Force displacement data correlation. 


A closer examination of the Case 3 model behavior is helpful to understand its 
functionality in more detail. A version of the Case 3 model for each mesh size was 
created where all elements up to 15.0 mm away from the delamination front were split 
and tied together (i.e., aw = 15.0 mm). To suit this condition, 15, 6, and 3 elements are 
split ahead of the delamination front in the 1.0 mm, 2.5 mm, and 5.0 mm mesh models, 
respectively. Using these modified Case 3 models, the force, F z , in each of the nodal ties 
in the split elements is normalized by the value at the crack tip and plotted in Figure 7b 
versus a. The definitions for a and cw are shown again in Figure 7a for convenience. 
Figure 7b shows that the applied DCB load is not carried entirely by one or two ties in 
the vicinity of the crack tip, as Case 2 effectively assumes, but by as many nodes as are 
present in a fixed distance ahead of the delamination front. 

The tie force in all mesh sizes is seen in Figure 7b to converge to zero at 
approximately a = 10.0 mm indicating that in this particular DCB specimen, all nodes 
within the range of 0 < a < 10.0 mm help to resist the applied load. Intuitively then, if 
there are not ties present within the 10.0 mm range or up until the force has sufficiently 
vanished, the applied load is resisted globally by what amounts to an inaccurate 
boundary condition for the DCB bending arms. This is what happens in Case 2 where 
there is only one node resisting the applied load ahead of the delamination. Figure 7b 
also illustrates why the 5.0 mm mesh is the most accurate within Case 2 as the single 
supporting node ahead of the delamination is located further along the curve closer to 
the point where the tie loads vanish. At the time of this paper, a generally applicable 
method to determine the distance required to split elements ahead of the delamination 
has not been determined. 


Figure 7c compliments the previous conclusions made using Figure 7b. Energy 
release rate convergence for each Case 3 mesh is plotted versus a m ax. Note that the first 
data point for each curve in this plot is equivalent to Case 2. In all curves, energy release 
rate is obtained from a node halfway across the specimen width at a common, but 
arbitrary, load step prior to any damage initiation. Energy release rate predicted by all 
meshes converges to the same value which also matches an analytical solution [24] 
given by 


G,= 


P 2 a 2 

BEI 


(9) 


where P is the applied load, a is the initial crack length, B is the specimen width, E is 
the stiffness modulus set equal to En here, and / is the area moment of inertia. 
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Figure 7. Tie forces in split elements ahead of the crack tip and Gi as a 
function of distance elements are split ahead of delamination. 


The remaining results are shown for Case 3 (as described in Table I) only. Figure 8 
is a plot of energy release rate normalized by G c across the delamination front. The data 
shown correspond to the load increment just before the first instance of damage growth. 
Results from each mesh size are included along with an analytical solution given by 
equation (9). While some small variance can be seen, there does not appear to be an 
appreciable amount of mesh dependence on energy release rate calculation for the mesh 
sizes considered. An exception to this is at the edges where the resolution in the finer 
meshes allows the drop in G to be characterized in more detail. The average Gi from the 
1.0 mm mesh model has an error of approximately 5% when compared to the analytical 
solution. This amount of error agrees with a similar observation made in the force 


displacement curves in Figure 6c. Based on Figure 6c and Figure 8, meshes using the 
Case 3 conditions sized corresponding to a Aa/h ratio range of 0.66 - 3.33, where Aa is 
the mesh size and h is half of the total laminate thickness, will produce accurate energy 
release rates in the enriched element for Mode I delaminations. 
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Figure 8. DCB: Results from mesh convergence study for normali z ed energy release 
rate, Gt/Gi c , distribution across delamination front (a = a 0 = 30.5 mm). 


In Figure 9a, the DCB model with a 1.0 mm mesh is shown in plan view at the end 
of a load increment identified as “Frame 1” just prior to the load step where Gt exceeds 
Gc for the first time. For illustrative purposes, a two-tone contour plot of deflection in 
the z-direction (i.e. out-of-plane) is shown where the lighter gray represents deflections 
greater than z = 0 (z = 0 lies on the specimen mid-plane). In this manner, elements that 
are delaminated can be identified by the fact that they have split and are deflected 
upwards in the positive z-direction. Following Frame 1, sequential load increments. 
Frames 2-5 seen in Figures 9b-9e, demonstrate the process for the delamination to 
advance one row of elements with ties at nodes across the delamination front being 
released according to their energy release rate exceeding G c as shown in the plots of Gt 
for each frame. This process repeats one row at a time throughout propagation. This 
manner of delamination advancement through the mesh is not exactly physically 
accurate, as typically delaminations in a DCB specimen form a curved front that 
advances at the same rate across the width of the specimen, however, it is equivalent to 
the behavior seen in the high-fidelity benchmark model [23], The large spikes in energy 
release rate located at comer nodes on the delamination front do not appear to cause any 
inaccuracies in the simulation of this specimen. 


See b) - e) 
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Figure 9. DCB: Delamination front geometry and energy release rate at the first 
instance of crack growth (a = a 0 = 30.5 mm). 


CONCLUSIONS 


A formulation was presented for a novel enriched shell finite element capable of 
delamination simulation in composite laminates. VCCT is combined with the Floating 
Node Method to simulate delamination growth at an arbitrary ply interface in a laminate 
layup. The element uses an adaptive splitting approach for damage characterization that 
allows for straightforward low-fidelity model creation and a numerically efficient 
solution. 

The enriched shell element was verified for Mode I delamination simulation by 
comparing results against a numerical benchmark study. The mesh configuration in the 
vicinity of the delamination front was investigated and it was determined that ties 
connecting split elements should not constrain rotations and that more than one element 
ahead of the delamination must be split and tied together. Following these guidelines, 
the enriched element was shown to exhibit little mesh dependency while accurately 
predicting structural stiffness and energy release rate. More specifically, meshes sized 
to a crack extension length-to-laminate thickness ratio, Aa/h, of roughly 0.66 - 3.33 
were all shown to result in accurate behavior. Investigation is ongoing in terms of 
mi xed-mode delamination simulation and extension of the element enrichment to 
include transverse/out-of-plane type damage features such as matrix cracks. 
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